arXiv:1312.0600vl [astro-ph.HE] 2 Dec 2013 


Accretion disks around binary black holes of unequal mass: 
GRMHD simulations near decoupling 

Roman Gold 1 , Vasileios Paschalidis 1 , Zachariah B. Etienne 1,2,3,4 , Stuart L. Shapiro 1,5 
1 Department of Physics, University of Illinois at Urbana- Champaign, Urbana, IL 61801 
2 Physics Department & Joint Space- Science Institute, 

University of Maryland, College Park, MD 207 f 2 
3 Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771 
4 Department of Mathematics, West Virginia University, Morgantown, WV 26506 
5 Department of Astronomy & NCSA, University of Illinois at Urbana- Champaign, Urbana, IL 61801 

and 

Harald P. Pfeiffer 6,7 

6 Canadian Institute for Theoretical Astrophysics, University of Toronto, ON M5S 3H8, Canada 
7 Canadian Institute for Advanced Research, Toronto, ON M5G 1Z8, Canada 

We report on simulations in general relativity of magnetized disks onto black hole binaries. We 
vary the binary mass ratio from 1:1 to 1:10 and evolve the systems when they orbit near the binary- 
disk decoupling radius. We compare (surface) density profiles, accretion rates (relative to a single, 
non-spinning black hole) , variability, effective o-stress levels and luminosities as functions of the mass 
ratio. We treat the disks in two limiting regimes: rapid radiative cooling and no radiative cooling. 
The magnetic field lines clearly reveal jets emerging from both black hole horizons and merging into 
one common jet at large distances. The magnetic fields give rise to much stronger shock heating 
than the pure hydrodynamic flows, completely alter the disk structure, and boost accretion rates 
and luminosities. Accretion streams near the horizons are among the densest structures; in fact, the 
1:10 no-cooling evolution results in a refilling of the cavity. The typical effective temperature in the 
bulk of the disk is ~ 10 5 (M/10 8 Mo) _1//4 (L/L e dd) 1 ^ 4 K yielding characteristic thermal frequencies 
~ 1O 15 (M/1O 8 M0) -1 / 4 (L/L e dd) 1//4 (1 + z) -1 Hz. These systems are thus promising targets for many 
extragalactic optical surveys, such as LSST, WFIRST, and PanSTARRS. 
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I. INTRODUCTION AND MOTIVATION 

Supermassive black hole (SMBH) binaries can form 
in magnetized plasma following galaxy mergers, via bar- 
mode instability in rapidly rotating supermassive stars, 
or by other dynamical processes [1]. After formation, a 
combination of dynamical friction and gas-driven migra- 
tion is likely to catalyze the binary inspiral into the grav- 
itational radiation-driven regime [1-7]. The exact details 
of these processes, including the “last-parsec problem”, 
remain active areas of research (see e.g. [8-10] and refer- 
ences therein). As a result, event rates and population 
synthesis studies at this stage are highly uncertain [11]. 

The exciting prospect of a simultaneous observation 
of both electromagnetic (EM) and gravitational waves 
(GWs) arising from accreting binary BHBHs makes these 
systems prime targets in the era of multi-messenger as- 
tronomy. Such observations will enable us to determine 
the binary masses, BH spins, redshift and even determine 
the Hubble constant to better than 1% [12-14]. 

Gravitational waves (GWs) from SMBHs are expected 
to be detected by planned GW interferometers such as 
eLISA/NGO detectors [15], sensitive to GW frequencies 
10 -5 — 1Hz, and the currently operating Pulsar Timing 
Arrays [16-18], sensitive to frequencies 10 -9 — 10 -6 Hz. 
As SMBHs are typically believed to have masses in the 
range 10 6 — 1O 9 M 0 , the GW strain at orbital separa- 
tion d = 10 M - the value adopted in this work - is 


h ~ lO- 16 (d/lOM)- 1 (M/lO 8 M 0 )(D/16Gpc)- 1 . Here 
M is the total mass of the binary, and we normal- 
ized to a luminosity distance D corresponding to red- 
shift z ~ 2 in a standard ACDM Universe. The 
corresponding gravitational-wave frequency is /gw ^ 
10“ 4 (M/10 8 Mq) _1 (d/10M)~ 3/2 [(1 + z)/3]~ l Ez. The 
expected GW strain is above the eLISA strain sensitiv- 
ity at these frequencies [15], hence these systems will be 
detectable by eLISA. In particular, for M 10 6 m 0 , 
the value of /gw at d — 10M falls well within the 
eLISA sensitivity band even for larger redshifts, while 
for M ~ 10 8 Mq and large redshifts (z ~ 10), the value 
of /gw at d = 10 M is marginally within the eLISA sen- 
sitivity band. However, the inspiral time from these sep- 
arations is tow ~ 2O(M/1O 8 M 0 ) days assuming equal 
masses. Hence, these systems will quickly enter the 
eLISA sweet spot, and EM precursor signals can trig- 
ger targeted GW searches with a convenient lead time of 
several days. 

While awaiting the first detection of GWs, currently 
operating and future electromagnetic (EM) detectors 
such as LSST [19], WFIRST [20] and PanSTARRS [21] 
are promising instruments to identify accreting BHBH 
systems in the EM spectrum. Important steps have al- 
ready been made toward realizing this goal. 

Currently, we know of one spatially resolved SMBH 
binary candidate at an orbital separation d ~ 7pc: 
0402+379 [22]. Other spatially unresolved, SMBH bi- 


2 


nary candidates have been found, including OJ287 [23- 
26] and SDSS J1536+0441 [27, 28]. Binary AGN can- 
didates have been singled out based on offsets in the 
broad line and narrow line regions, emission line pro- 
files, and time variability [29, 30]. Recently, very- long 
baseline interferometric observations interpreted ejection 
components from AGN cores as undulations caused by 
the precession of the accretion disks around a SMBH bi- 
nary [31]. A simplified model was applied to two AGN 
sources; for 1823+568 their analysis yields d ~ 0.42pc, 
and a mass ratio q in the range 0.095 < q < 0.25, while 
for 3C 279 d ~ 2.7pc, q ~ 0.36. However, given the 
lack of a robust circumbinary accretion disk theory these 
results are at best preliminary. Nevertheless, they still 
motivate a study of accretion flows onto supermassive 
BH binaries with different mass ratios, which we initiate 
in this work in general relativity (GR). 

Other features identified as “smoking guns” for binary 
black holes include BH recoil/kicks [32], used to explain 
the large velocity offsets between emission lines in AGN 
spectra, as well as observed kinks in jets probably due 
to changes in BH spin (X-shaped radio sources), a past 
merger event, or precession effects [33-35]. Modifications 
to the line profiles have also been proposed as promis- 
ing characteristic features to distinguish binary BH AGN 
from classical, single BH AGN sources [36]. 

To assist and solidify all these detection efforts, it is 
crucial to identify and model possible electromagnetic 
“precursor” and “afterglow” signatures [37-43]. 

Depending on the physical regime the properties char- 
acterizing the gas can differ considerably, and different 
accretion models are applicable. For example, if the gas 
has little angular momentum, the accretion flow resem- 
bles the binary analog of a Bondi-Hoyle-Lyttleton solu- 
tion [44-46] (see [47-49] for GR studies). If the gas has 
significant angular momentum, then the gas can become 
rotationally supported and form a disk. 

For a BHBH embedded in a disk, one can identify sev- 
eral different regimes based on the time scales for mi- 
gration of the binary. For a SMBH binary engulfed by 
a (thin) disk at large separation, the migration of the 
binary is initially governed by binary-disk angular mo- 
mentum exchange mediated by (effective) viscosity [3]. 
At large enough separations, the reduced mass fi of the 
binary is less than the local interior disk mass (47r d 2 X 
where X is the surface density of the disk). This leads 
to the so-called disk-dominated type II migration occur- 
ring on the viscous time scale t v i s . As the migration 
proceeds, the reduced mass of the system becomes larger 
than the local disk mass and the migration enters the 
secondary-dominated type II regime, which occurs on a 
longer time scale fe a = (47r d 2 X//i) -/c x t v [ s > t v i S , where 
k a constant of order 0.4 if 47rd 2 X//i < 1 and 0 oth- 
erwise. Ultimately, the binary enters a regime at smaller 
orbital separations where angular momentum losses due 
to GWs dominate, and the binary migrates on the GW 
time scale £gw- regimes, the disk moves inwards 

on the viscous time scale. When t v [ s < fe d < tew or 


t V is < few < fed the disk can follow the inspiral of the 
binary and settle in a quasi-steady state - this is called 
the predecoupling regime. When £qw < feis? the binary 
decouples from the disk, i.e. the inward migration of the 
binary out-paces the inward drift of the disk. 

In this paper we focus on the phase near decoupling , 
while the post decoupling regime will be the subject of a 
future paper. We note that unlike eLISA, Pulsar Timing 
Arrays are sensitive only to SMBH binaries well within 
the predecoupling regime. 

Accretion onto a single BH has been studied in great 
detail for decades, and magnetohydrodynamic studies in 
GR have drastically improved our understanding of these 
flows (see [50] for a recent review). Many different disk 
models have been proposed in the literature. These mod- 
els range from geometrically-thin, optically thick disks 
[51, 52] and slim disks [53], to geometrically thick, opti- 
cally thin, radiatively inefficient accretion flows (RIAF) 
[54-58]. However, our understanding of accretion flows 
onto BHBHs remains poor and studies of these systems 
are still in their infancy. 

The first analytic Newtonian model and smooth par- 
ticle hydrodynamic simulation of a circumbinary accre- 
tion disk was given in [59]. Since then, other Newto- 
nian (semi-) analytic studies [40-42, 60, 61] and hydrody- 
namic simulations in 2D [62-65], and 3D [66-70] have fol- 
lowed. Newtonian magnetohydro dynamic (MHD) simu- 
lations were presented in [71], and Post-Newtonian MHD 
simulations in [72]. Many of these earlier studies ex- 
cluded the binary and most of the inner cavity from 
the computational domain, introducing an artificial in- 
ner boundary condition. The importance of treating the 
inner regions self-consistently has been discussed in [69], 
and only full GR calculations can achieve this goal reli- 
ably. 

A “GR- hybrid” orbit- averaged model for thin disks, in 
which the viscous part is handled in GR and the tidal 
torques in Newtonian gravity was introduced in [43], and 
GR hydro dynamical simulations of accretion onto BHBH 
binaries - taking into account the dynamical spacetime - 
have been performed in [47, 73-75]. To date the only GR 
magnetohydrodynamic (GRMHD) simulations of disk ac- 
cretion onto BHBHs that account both for the dynamical 
spacetime and the BH horizons were presented in [76]. 

Using a different approach by assuming that the B-field 
is anchored to a circumbinary disk outside the computa- 
tional domain, [77-80] modeled EM signatures by solving 
the GR force- free electrodynamic equations. 

Close to merger a single BH remnant is formed on a 
time scale much shorter than the dynamical time scale 
in the disk. The mass and angular momentum of the 
remnant BH is different from the total mass and angular 
momentum prior to merger due to GW emission, causing 
a quasi-instantaneous perturbation to the disk. This ef- 
fect has been modeled using hydro dynamical and MHD 
simulations in Newtonian gravity and in GR [81-85]. 

A realistic and ideal 3D global model for a circumbi- 
nary disk around a SMBH binary near the decoupling 
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radius requires: a) a fully relativistic treatment of grav- 
itation in a dynamical spacetime, b) GRMHD for the 
plasma flow, c) realistic cooling processes, and d) radia- 
tive transfer in curved spacetime. Simulations incorpo- 
rating these effects must also have high resolution and 
long integration times (several viscous time scales). How- 
ever, including all these ingredients in one simulation 
would make these computations prohibitive, because the 
wall-clock times required to integrate for even ^ 5 vis- 
cous time scales at the inner disk radius at high resolution 
are far too long with current supercomputer resources. 
Thus, some of these ingredients must be relaxed in order 
to obtain a qualitative understanding of these systems. 
For this reason, the models in this work feature a) and 
b). High resolution is only adopted for a few models. In 
addition, we model radiative cooling by a simple cool- 
ing function and consider the extreme opposite limits of 
“rapid” cooling and “no cooling” to bracket the possibil- 
ities. 

In this paper we study the effects of the binary mass 
ratio on the disk near decoupling. We vary the BHBH 
binary mass ratio q = Mi/ M 2 < 1, considering 1:1, 1:2, 
1:4, 1:8 and 1:10 mass ratios. The mass ratio regime 
0-1 < g < 1 is shown to be of high astrophysical relevance 
in [86, 87]. Also, it has been argued that BHBHs forming 
in major galaxy mergers will typically have mass ratios 
in the range (0.01,1) [88]. These results motivate our 
choice of mass ratios. Note that Newtonian simulations 
studying the effects of mass ratio in the context of 2D 
thin-disk models were performed in [64, 65]. 

The new aspects of our work are the inclusion of rela- 
tivistic gravitation to resolve the crucial physics near the 
BH horizons, effective viscosity arising from the magne- 
torotational instability [89, 90], geometrically thick disks, 
three spatial dimensions (3D), a T-law equation of state 
(EOS), and effective cooling. We model the decoupling 
epoch by fixing the binary separation d, evolving the 
spacetime by rotating our conformal-thin-sandwich ini- 
tial data [91-94] as we have done in [76, 95]. The depen- 
dence of the flow variability, EM signatures, the magnetic 
field structure and the matter dynamics inside the low- 
density “hollow” on the mass ratio q are investigated. We 
estimate thermodynamic properties of the gas and scale 
our results with the binary total mass and the luminos- 
ity in units of the Eddington luminosity where feasible. 
From this, emission characteristics including typical ther- 
mal frequencies and luminosities are given and relevant 
detectors are discussed. 

In the absence of any observational constraints on the 
thermodynamic state of accreting BHBHs, the models we 
consider in this work adopt an adiabatic EOS governed 
by an adiabatic index T = 4/3, appropriate for radiation 
pressure-dominated, optically thick disks. This choice 
is motivated in part both by theory and observations. 
First, accretion disk theory [50] predicts that the inner 
part of circumbinary disks around SMBH binaries are op- 
tically thick, radiation pressure dominated for a large set 
of possible disk parameters (see e.g. [4, 52] and the next 


section). Furthermore, as discussed in [3], steady-state 
disk models predict that radiation pressure-dominated 
circumbinary disks will channel more material into the 
cavity for a given central mass. This finding makes radi- 
ation pressure-dominated circumbinary disks promising 
sources for electromagnetic counterparts. 

Second, AGN data from the Sloan Digital Sky Survey 
(SDSS) [96-98], the AGN and Galaxy Evolution Survey 
(AGES) [99], XMM-COSMOS [100] and others [101] cov- 
ering mostly type I AGNs and the local Universe to red- 
shift of 2 < 5, reveal Eddington ratio distributions in the 
range 0.01 < L/L^dd ^ 1 with the tendency that higher 
L / -7/Edd values occur at higher z (and possibly smaller 
central mass M). A comparison of accretion time scales 
with the age of the Universe suggests that earlier ac- 
cretion episodes are closer to the Eddington limit, sim- 
ilar to the recently discovered quasar at z ~ 7.1 with 
a central mass of ill ~ 2 x 10 9 M 0 [102] accreting at 
L/L^dd ~ 1.2 ± 0.5. These surveys therefore motivate 
studies of disks accreting near the Eddington luminosity, 
for which radiation pressure is important. Moreover, ra- 
diation pressure-dominated disks accreting near the Ed- 
dington limit are more likely to be detectable, even at 
large cosmological redshifts. 

This paper is structured as follows: In Sec. II we 
present a qualitative discussion of the range of param- 
eters and the associated physical regimes for which our 
simulations are appropriate. In Sec. Ill we describe our 
methods and techniques for evolving the spacetime, fluid, 
and magnetic fields, as well as our simple cooling pre- 
scription. We also present the different cases we consider 
and list our diagnostics for characterizing the accretion 
flow and EM signatures. In Sec. IV we show several tests 
we performed to motivate our numerical setup. In Sec. V 
we present the results from our simulations. We conclude 
in Sec. VI with a summary of the main results and a dis- 
cussion of future work. 

Here and throughout we adopt geometrized units, 
where G — c — 1, unless otherwise stated. 


II. QUALITATIVE CONSIDERATIONS 

In this section we use the Shakura-Sunyaev/Novikov- 
Thorne thin-disk model [51, 52] to make rough analytic 
estimates regarding the physical regime our disk mod- 
els probe. While the S hakur a- Sunyaev/ Novikov- Thorne 
model strictly applies to a viscous flow onto a single BH 
neglecting the binary tidal torques, we expect that our 
qualitative analysis will apply roughly to our circumbi- 
nary disk models, which have H/R = 0(0.1), where H 
is the disk scale height. This expectation should be best 
realized outside the binary orbit because the ratio of the 
tidal to viscous torques decays quickly far away from the 
binary orbit (see e.g. Fig. 6 in [43]). 
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A. Radiation pressure dominance 

1. Shakura-Sunyaev/ Novikov- Thorne model 

The simulations reported here apply to any total 
(ADM) binary black hole mass M and, since we neglect 
the self-gravity of the gas, to any rest-mass density po? 
provided it has the same initial disk profile adopted here. 
The quantities that are fixed, in addition to the initial 
disk profile, are the adiabatic index T appearing in the 
adopted ideal gas EOS, the ratio of the initial magnetic- 
to-total disk pressure, the initial magnetic field profile, 
the initial binary separation and the cooling law. In 
this section we use the steady-state thin-disk solution for 
a Shakura-Sunyaev /Novikov- Thorne disk about a single 
BH of mass M to estimate the physical values of some 
of the gas dynamical quantities as functions of M and 
luminosity L. We show below that for a range of astro- 
physically relevant choices of these parameters the disk 
is thermal radiation pressure-dominated, and this fact 
motivates our setting T = 4/3. 

Neglecting the perturbative role of the secondary tidal 
torque, the steady-state accretion flow in a geometri- 
cally thin disk is uniquely specified by the Shakura- 
Sunyaev viscosity parameter <a, the central mass M 
and the accretion rate M. The quantity M is spec- 
ified in turn by the disk luminosity L and efficiency 
5. The Shakura-Sunyaev/Novikov-Thorne model [51, 52] 
describes a disk that is radiation pressure-dominated in- 
side a radius ri nner . In this region the radiation to gas 
pressure ratio is [50] 



where e = L/ Me 2 is the radiative efficiency, L^dd — 1-3 x 
10 46 (M/10 8 M©)erg/s, and we chose the normalization 
for the a parameter based on typical values found in our 
simulations. For a thin-disk model the efficiency ranges 
from 0.057 for a non-spinning BH to 0.42 for a maximally 
spinning BH, so we scale to a value residing between these 
limits. The size of the inner region is determined by the 
condition ifed/^gas = 1, for which Eq. (1) yields 



The geometrically thick disks we evolve in this work 
extend radially out to r ou t ~ 100M— 200M. Hence, when 
scaling the accretion rate such that our models accrete 
near the Eddington limit, our models are fully immersed 
in this inner radiation pressure-dominated region. 


In this region the typical rest mass densities are [50] 


PO ~ 5.5 • 10" 12 



( M \f M \ g 
\1O 8 M0 ) y2.3MQ/yr J cm 3 


~ 5.5 • 10 -12 



m y 1 /_t_ y 2 /_e_y j© 

1O 8 M 0 ; UedJ V0.1/ cm 3 ' 


(3) 


(4) 


As we will show later, these characteristic values for 
Prad/^gas and po are comparable to the values found in 
our simulations. 

The Shakura-Sunyaev/Novikov-Thorne model predicts 
that the effective optical depth to absorption is r* ~ 0.02 
adopting the same normalizations as in the equation 
above. This well-known inconsistency near the Edding- 
ton limit of the Novikov- Thorne model is removed with 
the generalization to the slim disk model [53] . This model 
differs from a thin disk in that it allows for cooling to 
occur via advection, which dominates radiative cooling 
at high accretion rates (L > 0.3Z/Edd), thereby puffing 
up the disk. When scaling our models to accrete near 
the Eddington limit, they are closer to slim-disk models, 
which remain optically thick in this high luminosity limit 
[50, 53]. 

As we discuss later, near the Eddington limit and 
M = 1O 8 M 0 , the effective optical depth satisfies r* > 1 
in the bulk of our disk models, implying that the gas 
is optically thick to absorption and the photons eventu- 
ally thermalize. Thus, these qualitative considerations 
motivate the adoption of an adiabatic index T = 4/3, 
appropriate for thermal radiation pressure dominance. 

Note also that alternative disk models have been pro- 
posed, e.g. [4]. However, they largely share the prediction 
that the inner regions of the disk are radiation pressure 
dominated. 


2. Decoupling radius 


We estimate the decoupling radius ad by equating 
few = fes and solve for the separation to find [76] 




2/5 


(5) 


where we assumed that the inner disk edge radius set- 
tles to ri n ~ 1.5 d as typically found in our simulations, 
a is the Shakura-Sunyaev turbulent viscosity parame- 
ter, and fj = 4g/(l + q ) 2 (see also [5]). Notice that in 
contrast to geometrically thick disks, where the decou- 
pling radius is a few tens of M, geometrically thin disks 
have ad/M ~ 100. The decoupling radius estimate (5) 
for the mass ratios considered in this work ranges from 
ad{Q = 0.1) ~ 8.5M to ad(q = 1) ~ 13. 3M. The normal- 
izations in Eq. (5) are based on typical values obtained 
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from our simulations. We choose the binary separation 
d ~ 10 M for all cases studied in this work, a value which 
is consistent with the crude estimate of Eq. (5). In the 
future we intend to start our evolutions at larger separa- 
tions in order to dynamically determine the decoupling 
radius and evolve through it as in [43]. 


III. METHODS 

The models we adopt here assume: a) circular binary 
orbits, neglecting the binary inspiral (justified for large 
separations; see Sec. II A 2), b) non self-gravitating disks, 
which likely is a good assumption (see, e.g. [103]), c) ideal 
MHD, d) no radiative feedback, e) an effective cooling 
scheme that brackets no cooling and rapid cooling. 

Some of these assumptions may not be obeyed strictly, 
e.g. the binary may become eccentric in the predecou- 
pling regime [6, 66-69, 104], or radiative feedback may 
become important near Eddington accretion rates. How- 
ever, our simulations still provide a qualitative under- 
standing of the physics that will be useful in designing 
the next set of more realistic models of binary black holes 
immersed in circumbinary disks. In this section, we de- 
scribe the initial data and computational methods we 
adopt to account for a)-e). 


A. Initial data and AMR hierarchy 

1. Metric initial data 

At large separation the binary inspiral time scale is 
much longer than the binary orbital period and the vis- 
cous time scale at the inner edge of the disk just be- 
yond the binary orbit. Accordingly, the inspiral can 
be neglected over many orbital periods. To model this 
epoch in GR, we adopt quasi-equilibrium conformal-thin- 
sandwich (CTS) solutions for the black hole binary [91- 
94] . The spacetime initial data satisfying the CTS equa- 
tions correspond to a circular orbit and possess a helical 
Killing vector. The CTS initial data have been generated 
using the spectral techniques described in [105] as imple- 
mented in the Spectral Einstein Code (SpEC) [106] (see 
also [107]). We list the initial data parameters describing 
our spacetimes in Table I. 


2. Matter and B-field initial data 

For the fluid we use the same family of equilibrium 

disk models around single BHs as in [95, 108, 109]. We 
choose the initial inner disk edge radius r^ n? o = 18M and 
specific angular momentum /(r^o) = 5.15 M 2 * around a 
non-spinning BH as in [76]. However, the disk model is 


TABLE I. CTS initial data parameters for the BHBH vac- 
uum spacetime. Columns show mass ratio ( q ), ADM mass 
(Madm) and angular momentum (Jadm), and irreducible 
masses ( M - rr ), and apparent horizon radii ( r l hor ) for the two 
black holes. Diagnostics generating these quantities, but com- 
puted from vacuum, test simulations agree with these values 
to within one part in 10 4 . 

q Madm Iadm Mf TT M^ r r h or r^ or 

1 : 1 0.98989 0.96865 0.50000 0.50000 0.42958 0.42958 

1 : 2 0.99097 0.85933 0.66667 0.33333 0.60192 0.27312 

1 : 4 0.99342 0.61603 0.80000 0.20000 0.75140 0.15832 

1 : 8 0.99589 0.37868 0.88889 0.11111 0.85640 0.08618 

1 : 10 0.99656 0.31652 0.90909 0.09091 0.88081 0.07022 


not identical to [76] due to the different poly tropic index 
(T = 4/3 here versus T = 5/3 in [76]). 

We seed the initial disk with a small, purely poloidal B- 
field using the same procedure as in [76, 110]; see Fig. 1. 
The field is dynamically unimportant initially: the initial 
maximum value for the ratio of magnetic Pm to total 
pressure P is 0.025. All cases we consider in this work 
are initialized with the same disk and magnetic field. 



FIG. 1. B-field lines (white curves) and volume rendering of 
rest-mass density normalized to its maximum value at t — 0. 
The two black dots at the center indicate the BH apparent 
horizons for the q— 1 case. 

Although we will qualitatively discuss how the evolu- 
tion of these matter initial data depends on the binary 
mass ratio, it is important to stress that the goal of our 
work is to assess how the final relaxed state of the disk de- 
pends on the binary mass ratio. The initial data, which 
governs the early evolution, has no physical significance. 

B. Evolution equations and techniques 

1. GRMHD evolution 

We use the Illinois GRMHD adaptive-mesh-refinement 
(AMR) code [111-113], which adopts the Cactus/Carpet 
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infrastructure [114-116]. This code has been extensively 
tested and used in the past to study numerous sys- 
tems involving compact objects and/or magnetic fields 
(see e.g. [110, 117-121]), including black hole binaries in 
gaseous media [47, 76, 95]. 

The code solves the equations of ideal GRMHD in a 
flux-conservative formulation [see Eqs. (27)-(29) in [112]] 
employing a high-resolution-shock-capturing scheme (see 
[111, 112] for details), and including effective cooling 
source terms [see Eqs. (65) and (66) in [117]]. To en- 
force the zero- divergence constraint of the magnetic field, 
we solve the magnetic induction equation using a vec- 
tor potential formulation [see Eq. (9) in [113]]. As our 
EM gauge choice we use the generalized Lorenz gauge 
condition we developed in [76] and used in [118, 121]. 
We choose the damping parameter £ = 8/M. The 
advantage of this gauge condition is that it leads to 
damped , propagating EM gauge waves preventing spu- 
rious magnetic fields from arising near AMR boundaries 
even more effectively than the standard Lorenz gauge 
choice (£ = 0) [113]. The damping properties of the gen- 
eralized Lorenz gauge are crucial for stable and accurate 
long-term GRMHD evolutions. The “algebraic” gauge 
condition used in the first GRMHD simulations adopt- 
ing A-field evolution (see e.g. [49, 112, 122]) was shown 
in [113] to suffer from spurious conversion of EM gauge 
modes into physical modes and vice-versa, due to inter- 
polation at AMR boundaries. These spurious magnetic 
fields contaminate the evolution and the effect is exacer- 
bated when matter crosses refinement boundaries. 

To close the system of equations we use a T-law EOS 

P = p 0 e(T - 1), (6) 

where P is the total pressure, po the rest-mass density, 
and e the specific internal energy. This EOS allows for 
shock heating. We choose T = 4/3 appropriate for radi- 
ation pressure-dominated, optically thick disks. 


2. Metric evolution 

The metric evolution is treated under the approxima- 
tion that the inspiral time scale due to GW emission is 
long compared to both the binary orbital period and the 
viscous time scale of the disk. Hence, we can neglect 
the inspiral for multiple binary orbits. The CTS initial 
data we adopt possess a helical Killing vector, which im- 
plies that the gravitational fields are stationary in a frame 
corotating with the binary. As a result, we can perform 
the metric evolution in the center-of-mass frame of the 
binary by simply rotating the initial spacetime fields as 
was done in [47, 76]. This technique simplifies our com- 
putations substantially. In addition, the rotating metric 
method facilitates our evolving dynamically to relaxed 
disk/magnetic field initial data for the inspiral. 

To implement this method, we map the CTS solution 
from the spectral grid onto three grids corresponding to 
three partially overlapping spherical coordinate systems: 


One spherical coordinate system covers the entire evolu- 
tion domain and is centered on the binary center of mass, 
and two smaller ones are centered on each BH. These new 
grids employ a logarithmic radial coordinate. We use the 
CTS solution stored on these spherical grids to interpo- 
late the data onto the Cartesian evolution grids whenever 
we perform the rotation transformation. 

We have checked that the mapping from the spectral 
grid to the spherical grids is implemented correctly by 
performing vacuum simulations that use the CTS solu- 
tion stored in the spherical grids as initial data. More 
specifically, we have computed several diagnostic quan- 
tities which characterize the BHs and the global space- 
time and compared them with the values known from 
the spectral CTS initial data (see Table I). These agree 
to within 1 part in 10 4 . Moreover, we have verified that 
a crude estimate for the orbital frequency of the binary 
(orbital trajectory traverses a full phase of 2tt) as deter- 
mined by a dynamical vacuum evolution agrees with the 
value given by the initial data (M9 = 0.028) to within 
^ 10%. We have computed the normalized L 2 ,n norm 
of the Hamiltonian and momentum constraint violations 
as introduced in Eqs. (59) and (60) in [123], with the 
modification that we split up the Laplacian operator into 
its individual components when computing normalized 
norms. We find the normalized norm of the Hamiltonian 
constraint to be dominant, with L 2 ,n(J~L) ~ 2%. We con- 
clude that the CTS solutions are mapped correctly and 
accurately. 


3. Cooling 


Without cooling, the binary tidal and the viscous 
torques act to gradually heat and puff up the disk. “Ad- 
vective” cooling, which is crucial in slim-disk and ADAF 
models [50] , is self-consistently accounted for in our simu- 
lations. However, adding radiative cooling may be neces- 
sary to achieve a steady state. Realistic radiative cooling 
based on actual physical mechanisms depends on com- 
plicated microphysics, which we do not model here, but 
intend to incorporate in future studies. To model steady- 
state solutions in this work, we introduce “radiative” 
cooling via an effective cooling leakage scheme. This 
scheme is, strictly speaking, valid in the optically thin 
regime. While this is a very crude approach to radiative 
cooling, treating both extremely rapid radiative cooling 
as well as no radiative cooling can help bracket the pos- 
sibilities. Different formulae for the cooling emissivity A 
have been proposed in the literature: 

I. Non-exponential cooling [72, 124]: 


A = 


Pot /A S + 

Tcoo\ So 



(7) 


where S = K = is an entropy parameter, and 
P 0 

S 0 the initial target value. We call this emissiv- 
ity “non-exponential” , because the effective cooling 
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time scale for this scheme is not just r coo i, but de- 
pends on the internal energy e (see Appendix A). 

II. Exponential cooling [117, 119] 

\ Po^th Po^o {AS | AS jr ^ 

A - “ FZ K~q~ + \ W 

T cool ^Tcool ^0 ^0 

where 6$ is the internal energy calculated using So, 
and eth = e — eo- In this scheme the effective cooling 
time scale is r coo i- 

Both emissivities dissipate all shock-induced thermal en- 
ergy, driving the entropy of the gas to its initial value. 
We use prescription II, instead of I, because we have 
found that prescription I is prone to the development of 
a Courant instability, as the effective cooling time scale 
of this scheme depends on the amount of shock heating, 
which can be very strong in low-density regions (see e.g. 
Appendix B in [125]). Thus, to stabilize the simulations 
with prescription I, one typically excludes cooling of the 
low-density regions or unbound matter [72, 76]. As both 
the BH horizons and the low-density cavity is included 
in our computational domain (unlike earlier studies), we 
find that the strong shock heating inside the cavity in 
conjunction with emissivity I leads to a numerical insta- 
bility. In order to bracket the effect of cooling, this inner 
cavity needs to be cooled when cooling is enabled. In 
Appendix A, we present an analytic calculation illustrat- 
ing the above considerations. We demonstrate that shock 
heating of matter in the cavity is important in Sec. IV C. 

To model “rapid” radiative cooling we set the cooling 
time scale equal to 10% of the local, Keplerian time scale 
t coo \/M = 0. It K ep / Af = 0.1 • 27r(r/M) 3 / 2 , where r is the 
cylindrical radial coordinate measured from the center of 
mass of the binary. In order to prevent the cooling time 
from becoming prohibitively small as r 0 we floor the 
cooling time at r coo i > 10M. Throughout this paper we 
refer to cases with A % 0 as the cooling cases and A = 0 
as the no-cooling cases. 


4- Evolution grids & models 

We use a hierarchical, box-in-box adaptive mesh pro- 
vided by the Cactus/Carpet infrastructure [115, 116]. We 

constructed two sets of nested boxes, with one set cen- 
tered on each BH, on which we discretize the GRMHD 
evolution equations. The coarsest level has an outer 
boundary at r = 240M. Due to a range of resolution 
requirements related to the different sizes of the BHs, 
different models use different number of refinement lev- 
els, which in turn yields different finest grid spacings (see 
Table II). We set up the locations of our AMR bound- 
aries such that the computational grids resolve both the 
BHs and the inner disk region for the given resources. 


The grid spacing is also motivated by both MRI resolu- 
tion requirements and the results gleaned from test runs 
involving hydrodynamic disk evolutions around a single, 
non-spinning BH, which we report on in Sec. IV. 

In Table II we also list the distinguishing characteris- 
tics of the different cases we consider in this work. The 
labels are chosen to designate the mass ratio, whether 
cooling is applied or not, and the resolution. For ex- 
ample, the label l:lnc-hr means mass ratio q = 1, no- 
cooling, and high resolution. 

We point out that the disparity in length scales (hori- 
zon vs. disk size) and time scales (the Courant condi- 
tion vs. viscous time scale) intrinsic to the circumbinary 
BHs disk problem introduces a large computational cost. 
Most of our simulations were run continuously (excluding 
queue waiting times) for ^ 2 months on Blue Waters, 
Kraken, Lone star, as well as the Illinois Relativity group 
36-node Beowulf cluster. The CPU hours used depended 
on the computer cluster, but we estimate that the simu- 
lations presented here required ^2x 10 4 * 6 CPU hours. 


C. Diagnostics 

We distinguish two types of diagnostics. The first type 
consists of probes of the MHD flow, including the den- 
sity and velocity profiles, accretion rates, luminosity esti- 
mates, magnetic field profiles, the establishment of MRI 
turbulence and jets, etc. The second type concerns prop- 
erties of the plasma such as local temperatures, optical 
depths, characteristic frequencies of emitted radiation, 
etc. The first type are straightforward to calculate from 
our simulation data, as they depend on the overall MHD 
behavior of the disk, and once we have chosen an EOS 
and cooling prescription, are independent of detailed mi- 
crophysics. We are quite careful and confident in re- 
porting these diagnostic quantities. The second type de- 
pends crucially on the specific physical values we assign 
to our nondimensional input parameters (e.g. BH mass 
and disk density) and to the microphysics that is not 
incorporated in our calculation (e.g. realistic radiative 
cooling and radiative transport). Nevertheless, we can 
make crude estimates for the latter quantities, and do so 
below. Although considerable caution must be applied to 
these estimates, they may serve as useful guides to subse- 
quent, more detailed investigations and to astronomical 
instruments that may be able to observe the scenarios we 
are simulating. 

The first type of diagnostics includes: 1) Accretion 
rate M as defined in [47]. We compute the total accre- 
tion rate onto the binary, and also the accretion rate onto 
each individual BH. 2) Fourier analysis of the accretion 
rate FFT(M ), targeted to identify possible quasiperi- 
odic signatures in the accretion flow. 3) Surface den- 
sity profile E(r) as defined in [95]. This diagnostic is 
also useful to compare with studies of ID orbit-averaged 
disk models. 4) Shakura-Sunyaev ^-stress parameter 


TABLE II. List of grid parameters for all models. Equatorial symmetry is imposed in all cases. The computational mesh 
consists of two sets of nested AMR grids, one centered on each BH, with the outer boundary at 240M in all cases. From left 
to right the columns indicate the case label, mass ratio g, whether cooling is included or not, the coarsest grid spacing Ax max , 
number of AMR levels around the primary (BH) and the secondary (bh), and the half length of each AMR box centered on 
each BH. The grid spacing of all other levels is Ax max /2 n_1 , n — 1,2,..., n max , where n is the level number such that n — 1 
corresponds to the coarsest level. A dash indicates “no information available”. 


Case 

Q 

cooling? A^ max 

levels(BH) 

levels (bh) 


Grid hierarchy 

l:lnc-hr 

1:1 

No 

4.8 M 

7 

7 

240M/2"- 1 , 

n — 2, . 

..5, 240M/2 n , n = 6,7 

l:lnc-mr 

1:1 

No 

6.0 M 

7 

7 

240M/2 n_1 , 

n — 2, . 

..5, 240M/2 n , n = 6, 7 

l:lnc-lr 

1:1 

No 

8.0 M 

7 

7 

240M/2 n_1 , 

n — 2, . 

..5, 240M/2 n , n = 6,7 

l:2nc-lr 

1:2 

No 

8.0 M 

7 

8 

240M/2 n_1 , 

n — 2, . 

. .5, 240M/2 n , n = 6,7,8 

l:4nc-lr 

1:4 

No 

8.0 M 

7 

9 

240M/2 n_1 , 

n — 2, . 

..5, 240M/2 n , n = 6, . . . , 9 

l:8nc-lr 

1:8 

No 

8.0 M 

6 

10 

240M/2 n_1 , 

n — 2, . 

..5, 240M/2 n , n = 6,..., 10 

l:10nc-lr 

1:10 

No 

8.0 M 

6 

10 

240M/2 71-1 , 

n — 2, . 

..5, 240M/2 71 , n = 6,..., 10 

Onc-hr 

0 

No 

4.8 M 

6 

- 

240M/2 71-1 , 

n — 2, . 

..5, 240M/2 n , n = 6 

Onc-mr 

0 

No 

6.0 M 

6 

- 

240M/2 n_1 , 

n — 2, . 

..5, 240M/2 n , n = 6 

Onc-lr 

0 

No 

8.0 M 

6 

- 

240M/2 71-1 , 

n — 2, . 

..5, 240M/2 n , n = 6 

l:lc-mr 

1:1 

Yes 

6.0 M 

7 

7 

240M/2 n_1 , 

n — 2, . 

..5, 240M/2 n , 77 = 6,7 

l:2c-lr 

1:2 

Yes 

8.0 M 

7 

8 

240M/2 n_1 , 

n — 2, . 

. .5, 240M/2 n , n = 6,7,8 

l:4c-lr 

1:4 

Yes 

8.0 M 

7 

9 

240M/2 n_1 , 

n — 2, . 

..5, 240M/2 n , 77 = 6 ,..., 9 

l:8c-lr 

1:8 

Yes 

8.0 M 

6 

10 

240M/2 11-1 , 

n — 2, . 

..5, 240M/2 n , 7i = 6,..., 10 

l:10c-lr 

1:10 

Yes 

8.0 M 

6 

10 

240M/2 71-1 , 

n — 2, . 

..5, 240M/2 71 , n = 6,..., 10 

Oc-lr 

0 

No 

8.0 M 

6 

- 

240M/2 n_1 , 

n — 2, . 

..5, 240M/2 n , n = 6 


|jiEM ^ 

computed as a ~ a EM = ( ^ ) t where T™ is the 
dominant orthonormal component to the Maxwell stress- 
energy tensor evaluated using the tetrad defined in [126] 
(the brackets denote an orbit averaged quantity). More 
specifically, we report an azimuthally- and z- averaged 
a = a(r) profile, which can be used in ID orbit-averaged 
disk models. 5) Disk scale height H = Y,/po(z = 0). 
6 ) Inner disk edge radius r- in : In all E(r)-profiles we ob- 
serve that inside the cavity E(r) declines rapidly with 
decreasing r and becomes convex. We fit a fifth or- 
der polynomial to the orbit-averaged E(r) in the con- 
vex region at small r, and define r- in as the radius where 
the curvature of the E(r) fitting function is maximized, 
[E"(r)/( 1 + E^r) 2 ) 3 / 2 ]' = 0 where ' = d/dr. 7) EM 
Poynting luminosity (Lem) as defined in Eq. (1) of 
[121]. 8 ) Energy loss rate carried off by the outflowing 
gas L gas = § s To,/ gas ^dS. This surface integral must be 
performed in the asymptotically flat regime. Given that 
we do not perform the integration at an infinite radius, 
as a crude approximation to L gas we include in the in- 
tegration only matter that is unbound, i.e. , matter for 
which at large radii E = — uq > 1. We compute 7) and 
8 ) at several radii including 90M, 140M, 210M. 9) For 
cases where our cooling scheme is enabled, we compute 
the cooling luminosity L coo i = § s Tq ,( rad )<LS, which we 
estimate via the volume integral L coo i = f Auoa- s/ dyd 3 x. 
The volume integration is exactly equal to the surface in- 
tegration at steady-state and in spacetimes possessing a 
timelike Killing vector and when we ignore any radiation 
captured by the BHs. We also compute the bolometric 


radiative luminosity L 5 = Lem + L coo p 

The second type of diagnostics includes: 10) Opti- 
cal depth to true absorption r* = y/3r es T^ (Eddington 
approximation), where we assume pure hydrogen and 
where r es (tq) is the optical depth to electron scatter- 
ing (free-free absorption), calculated using the Thomp- 
son scattering opacity K es = 0.4cm 2 /g (Rosseland mean 
opacity % = 6.45 • 10 23 poL _ 3 - 5 cm 2 /g) as r es = K e S E 
(rff = ttffE) [127]. r* > 1 implies the matter is optically 
thick to absorption. 11 ) Local temperature of the matter, 
calculated by solving e = aT 4 /po + 2 &#TpoA%, where 
m p is the proton mass and the Boltzmann constant. 
12) In the cases with cooling, the effective disk temper- 
ature (in cooling cases), estimated by assuming that all 
cooling luminosity is emitted as black body radiation: 

T eff = [L cool /Ana{r 2 oal - rf n )] 1/4 , (9) 

where a is the Stefan-Boltzmann constant and r out is the 
disk outer radius. 13) Characteristic observed thermal 
radiation frequencies z/bb = /c^T e ff//i(l + z), where h is 
the Planck constant and z the cosmological redshift. This 
is calculated only when A ^ 0 . 14) Cyclotron emission. 
While we find the bulk of the disk to be optically thick 
near Eddington accretion rates, the low-density “cavity” 
is optically thin. From these regions cyclotron lines may 
be detectable. 15) In cases where A^O we compute the 
characteristic cyclotron frequencies v cy = eB/mc( 1 + z), 
where e is the electron charge, m the electron mass and 
B the magnitude of the magnetic field. 
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IV. TESTS AND RESOLUTION 
REQUIREMENTS 

In this section we describe tests we performed that 
motivate the choices for the grid resolution and cooling 
time scale. 


A. Hydrodynamical evolutions with B — 0 

To set a lower limit on the necessary resolution to 
perform our GRMHD simulations, we found the mini- 
mum resolution required so that our code maintains sta- 
ble equilibrium of an unmagnetized disk around a single 
non-spinning BH for several thousands of M of evolution. 
The equilibrium disk solution we use is described at the 
beginning of Sec. Ill A 2. Our study indicates that for 
the low resolution (Ax max = 8.0M) listed in Table II the 
surface density profile of the initial disk is maintained to 
within 2% throughout the disk for at least t ~ 5000M. 



—80 -60 -40 -20 0 20 40 60 80 


X/M 

FIG. 2. Contours of the AMRi-quality factor Q = Xyim/dx in 
the equatorial plane, corresponding to the equal mass (q = 1) 
medium resolution case [divide (multiply) by 1.33 (1.25) for 
the low (high) resolution case] at t — 0 M. We resolve the 
fastest growing MRI mode by > 10 gridpoints over a large 
part of the disk (the blue ring stems from extremely small 
values of Amri, when the vertical component of the B-field 
changes sign). 


B. MRI resolution requirements 

Here we check the conditions for MRI to operate in our 
disk models. For this to be the case three conditions have 
to be satisfied: (I) A magnetic field configuration must 
be present that violates the stability condition for MRI 
dVt/dr > 0, (II) The wavelength of the fastest-growing 
mode Amri has to be resolved by > 10 gridpoints [128- 
130], and (III) the B-field must be sufficiently weak for 
Amri ^ 2 H. In other words the wavelength of the fastest 
growing mode should fit in the disk [131]. 

Regarding (I) our initial disk model is unstable to the 
MRI because of the outwardly decreasing angular veloc- 
ity and the presence of an initially small poloidal mag- 
netic field. Regarding (II) we plot the quality factor 
Q = Amri /dx where dx is the local grid spacing (which 
jumps by a factor of two at AMR boundaries); see Fig. 2. 
One can see that we satisfy the criterion Q > 10 over a 
rather large portion of the disk initially except for the re- 
gion near the radius where the poloidal field changes sign 
and becomes very small. We chose our low-resolution 
grids such that condition (II) is satisfied at t = 0 in the 
innermost parts of the disk. 

Regarding (III) we plot a meridional (x, z)- slice of the 
rest-mass density overlayed by a line plot showing Amri 
as a function of x in Fig. 3. The plot shows that for 
the most part Amri fits within the disk. At the inner 
disk edge the MRI is likely to be suppressed initially, 
but as the evolution proceeds magnetic winding converts 
poloidal fields into toroidal ones lowering Amri and even- 
tually triggering the MRI near the inner disk-edge. 


i 
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0.6 
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0.4 

I 

Bo.o 

X/M 

FIG. 3. Rest-mass density contours (color coded) on a merid- 
ional slice, and Amri/2 (black solid line) at t — 0 M. The plot 
corresponds to equal mass (q = 1) but is the same among all 
cases considered in this work. For the most part Amri/2 fits 
within disk. 


C. Cooling 

We seek to compare two opposite limiting cases for 
each mass ratio: (I) No-cooling, r coo i t K ep for which 

A = 0. Here TRe P is the local Keplerian time scale which 
is comparable to the (shock) heating time scale; (II) ex- 
tremely rapid cooling, r coo i rx e p which we model with 
the effective emissivity A = po6th . 

^cool 
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To determine the value for r coo i at which cooling be- 
comes rapid (at least in the bulk of the disk), we per- 
formed the q = 1 cooling case using different cooling 
time scales. We concluded that rapid cooling requires 
a cooling time scale significantly shorter than the lo- 
cal Keplerian time scale. In Fig. 4, we plot K/Kq, 
where the entropy parameter S = K = -P/(po) an d 
So = Ko = K(t = 0) for a run with r coo i = O.lrKep 
(left panel), a run with r coo i = Tx ep (middle panel) and a 
run without cooling r coo i = oo (right panel). It becomes 
apparent that when r coo i = rx e p , K is not driven back to 
its initial value. Physically, this means that not all shock 
generated entropy is radiated away, hence r coo i = t K ep 
does not correspond to rapid cooling and steady state 
is not achieved. For r coo i = O.lrKep we find K/Kq ~ 1 
in the bulk of the disk. The values depart from unity 
only inside the cavity where low density gas exists and 
can be shock heated to very high K/Kq , demonstrating 
that shock heating is extremely strong in the low-density 
cavity. We adopt r coo i = O.lrKep m all our cooling simu- 
lations because it leads to rapid cooling, at least through- 
out the bulk of the disk. 


V. RESULTS 


In this section we present the results of our numerical 
simulations. First, in Sec. V A we show results from our 
resolution study. In Sec. VB, we directly compare the 
q = 1 binary case with B % 0 to the B = 0 case. Lastly, 
we report on the variation of our diagnostics with mass 
ratio for all magnetized cases in Sec. V C. 

Our simulations have two parameters that scale out 
of the problem: the total mass of the BHBH binary M 
and the rest-mass density of the disk. Alternatively, we 
can exchange one of these parameters with another pa- 
rameter that depends on these two free parameters. So, 
instead of the rest-mass density, in the results we quote 
we choose the Eddington ratio Tb/AEdd, where Lb is the 
bolometric EM luminosity described in Sec. IIIC. The 
relation between these parameters is determined as fol- 
lows: the accretion rate through the horizon must scale 
like M oc p 0 ,refAf 2 , where po,ref is a reference rest-mass 
density in the disk. We choose the maximum rest-mass 
density at t = 0 as the reference density, and our simula- 
tions determine the proportionality constant. For exam- 
ple, in the single, nonspinning BH case with cooling we 
find 


Pvei = 9 X 10 12 



M 

10 8 M© 


-1 

g cm -3 , 


where we have replaced the accretion rate via the follow- 
ing expression 


M = ibE" 1 = [ 


2.75 


V-^Edd 

U_\ ( M 


LftddS 


L-Rdd) V10 8 M g 


( 10 ) 


M© yr 


-1 


and where 5 = L^/M = 0.08 is the radiative efficiency 
as computed via our simulations for our adopted cooling 
law in the single, nonspinning BH case. In the no-cooling 
cases the only radiation luminosity estimate we have is 
the Poynting luminosity which is expected to be a small 
fraction of the total radiative luminosity. Hence, in the 
no-cooling cases we do not scale our results with the Ed- 
dington ratio. Instead, we choose a fiducial accretion rate 
similar to the one given in Eq. (10). 


A. Resolution study 

Here we present the results of our resolution study. 
For the single BH, no-cooling and BHBH equal mass, no- 
cooling cases we use the three resolutions (see Table II). 

In the single BH-case the average accretion rate varies 
little with resolution (see left panel in Fig. 5). The maxi- 
mum fractional difference of the mean accretion rate be- 
tween different resolutions is 15%. Other quantities show 
a similar behavior. These results indicate that the res- 
olutions used in this case are high enough for the main 
MRI effects to be captured and the results to be qualita- 
tively independent of resolution. However, the resolution 
is not sufficiently high to perform a formal convergence 
test. 

For the equal-mass case we observe a different behav- 
ior. The mean accretion rates appear to converge (see 
right panel in Fig. 5), but the low resolution run under- 
estimates the accretion rate by almost a factor of 2, while 
the medium and high resolution runs are in good agree- 
ment. For the latter resolutions mean accretion rates 
agree to within 9%. We conclude that for the q = 1 case 
our adopted medium and high resolutions are sufficient 
for drawing qualitative conclusions and that higher res- 
olutions are necessary for accurate quantitative results 
that reside in the convergent regime. 

The results of the resolution study in the q = 0 and 
q — 1 cases differ because of the distribution of matter in 
both cases and our grid setup. In the q = 0 case there is 
more matter close to the BH, where very high resolution 
refinement levels reside, whereas in the q — 1 case the 
bulk of the matter remains outside the inner edge of the 
disk, where the grid resolution is not as high. As we show 
later, this is not the case in our q < 1 models. There, 
more matter resides closer to the BHs, and hence closer 
to the high resolution levels. 

Based on our resolution study we conclude that the 
low resolution used in the equal-mass case is not suffi- 
ciently high to yield reliable results, but for the other 
mass ratios we consider in this work, our adopted low 
resolution suffices for a qualitative discussion. Thus, in 
the equal-mass cases results will be reported mainly from 
our medium resolution runs. 

We stress here that these simulations, which include 
all relativistic effects and resolve the BHs, are computa- 
tionally very demanding (see Sec. Ill) and increasing the 
resolution even further will incur a very large computa- 
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FIG. 4. K/Ko for three q — 1 runs with different r CO oi at t — 7700M. Left panel: r CO oi = O.lrKep- Middle panel: r CO oi = t K ep • 
Right panel: No-cooling (t coo i = oo). Cooling in the t coo i — TKep case cannot keep up with heating, causing K to drift 
away from its initial value. Values much smaller than TKep are required for rapid cooling. Steady state is not achieved in the 

Tcooi ^ "^Kepler Cases. 



FIG. 5. Left panel: rest-mass accretion rate vs time for the q — 0 no-cooling cases at low (lr), medium (mr), and high (hr) 
resolutions. Right panel: rest-mass accretion rate vs time for the q — 1 no-cooling cases at low, medium, and high resolutions. 
The accretion rate in the q = 0 (q = 1) case is normalized by the mean accretion of the high resolution q = 0 (q — 1) run. The 
horizontal lines indicate the mean accretion rate at low, medium and high resolution. 


tional cost. With increasing computer power and larger 
computer allocations we plan to improve our results in 
the near future. 


B. Significance of B-fields 

Previous hydrodynamic simulations and 
(semi) analytic models of circumbinary accretion 
disks using the simplified <a-disk model (e.g. [59, 63]) 
showed that the main feature in the equal-mass binary 
case is that the density inside and near the binary’s orbit 
remains substantially lower than in the single BH case 
(see, e.g., Fig. 6). Such an inner cavity or “hollow” can 
have important consequences for the emergent radiation, 
such as line emission due to small optical depth and 
small bolometric luminosity from the hollow. Any such 


difference between single BH vs. circumbinary accretion 
disks can provide a path to distinguishing a binary AGN 
versus a classical, standard (single BH) AGN [132]. The 
explanation for the existence of a hollow is that the 
binary tidal torques for q > 0.01, are strong enough to 
push most matter away from the binary orbit [3]. The 
effect is most prominent in geometrically thin disks, 
which arise when radiative cooling is highly efficient. 

However, even in the absence of viscosity or magnetic 
fields, the time-dependent tidal field strips off matter 
from the inner disk edge, giving rise to an accretion pat- 
tern consisting of two streams which penetrate the inner 
cavity and extend to the horizons of the BHs [3, 47, 73- 
75, 95]. 

Furthermore, recent MHD simulations [71, 72, 76] uni- 
versally revealed that the reduction of density inside the 
cavity in the binary case is not as substantial as previ- 
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ously thought. Such simulations explore regimes in which 
the disk is geometrically thick, which partially accounts 
for the difference. 

Here we present a comparison between magnetized and 
unmagnetized circumbinary accretion disks onto an equal 
mass BHBH, while all other physical and numerical pa- 
rameters remain identical, to illustrate the importance of 
magnetic fields in filling the hollow with dense material. 


1. Midplane- density 

Figure 6 demonstrates the striking differences between 
no-cooling evolutions with and without magnetic fields at 
the same integration time. Magnetic-free hydrodynamic 
evolutions severely underestimate both the density in the 
inner regions and the over density due to the spiral arms. 
This indicates that the amount of matter stripped off by 
tidal torques is small compared to the amount of matter 
flowing into the hollow due to MHD turbulence. 


2. T, -profiles, q — 1 

Next we compare the surface density profiles of mag- 
netized vs unmagnetized disks. We find different profiles 
between the two cases as shown in Fig. 7. The B = 0 
model remains relatively close to the initial data, apart 
from a slow, mild expansion due to tidal heating and 
shocks. When magnetic fields are present, the final disk 
profile is completely different from the initial data even 
though the binary torques are identical. This implies 
that the MRI-driven viscous torques have a much larger 
impact on the global disk structure than the binary tidal 
torques, except perhaps near the inner disk edge. 


3. Sensitivity to cooling 

We find a fundamental difference between B = 0 and 
B ^ 0 evolutions regarding their sensitivity to cooling. In 
the B = 0 case both A^O and A = 0 evolutions lead to 
essentially the same E-profile (see Fig. 7). This is in stark 
contrast to the 5/0 (magnetized) cases shown in Fig. 7, 
for which cooling has a strong impact, leading to a matter 
pile-up near the inner disk edge. As our particular choice 
of A serves to dissipate heat from shocks, we conclude 
that magnetic fields lead to far stronger shock heating in 
the disk than the binary tidal torques. 


Accretion rate 

The ratio of the time-averaged accretion rate with- 
out magnetic fields to that with magnetic fields is 
(Mbhbh,b=o) / (Mbhbh,b^o) < 1% (see also Sec. V C). 
This result applies to both cooling and no-cooling cases. 


In summary, B = 0 evolutions underestimate the ma- 
terial inside the cavity and accretion rates by orders 
of magnitude. Hence, incorporating magnetic fields is 
paramount for a proper treatment of circumbinary ac- 
cretion disks. 


C. Trend with mass ratio, 5/0 

In this section we discuss the dependence of our mul- 
tiple diagnostics on the binary mass ratio for our 5/0 
cases. We use results from the single non- spinning BH 
case (q = 0 ) to normalize and compare our results for the 
binary cases. 

1. Disk structure in the bulk and inside the cavity 

We begin this section by discussing the qualitative evo- 
lution of the disk rest-mass density po- For q 7 ^ 0 , the 
early evolution is similar for the different mass ratios, but 
departs strongly in the subsequent evolution depending 
on q. The onset of accretion occurs through two spiral 
streams, which remain attached to the horizons through- 
out the evolution, as shown in Figs. 8-15. Here we plot 
the rest-mass density contours in the equatorial plane. 
These streams are among the densest structures of the 
accretion flow, especially for the lowest non-zero mass 
ratio case q — 0.1 (see Figs. 14 and 15). Spiral density 
waves are launched near the inner edge of the disk, which 
propagate and dissipate into the outer disk. This feature 
can also be seen in Figs. 8-15 for all mass ratios. 

Late in the evolution we find that when q ^ 0, the 
density of the matter inside the “cavity” in A = 0 cases 
is larger than that in A / 0 cases (see Figs. 8 and 9). 
Hence, the amount of matter in the hollow is smaller 
when we allow rapid cooling. This dependence on cool- 
ing arises because of the larger disk thickness in the no- 
cooling cases, which leads to a reduction in tidal torques 
[5, 133] near the binary orbit, allowing matter to overflow 
more easily. 

We also find that the smaller the mass ratio the more 
matter pours into the cavity. This is anticipated from 
the Newtonian expression for the binary tidal torque [5, 
133], which decreases with decreasing q. As expected, the 
largest contrast arises between the q = 1 and q = 0 cases, 
which becomes obvious by simply inspecting the rest- 
mass density contours in the equatorial plane as shown 
in Fig. 16. One can see the main difference: The presence 
of a region of lowered density with two accretion streams 
near the BHs in the binary q = 1 case - the “hollow” (or 
“cavity” ) - and the absence of these features in the q — 0 
case (left panel). 

The relaxed disk structure in the predecoupling regime 
is not axisymmetric. The back-sloshing of material to- 
wards the inner disk edge occurs mainly in two opposing 
directions and leads to a gradual overdense feature in the 
disk, which has been referred to as a “lump” [71, 72], and 
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FIG. 6. Contours of rest-mass density normalized to the initial maximum rest-mass density in the equatorial plane at £ ~ 5000M 
for two q — 1 no-cooling cases. Left panel: B — 0. Right panel: B ^ 0. Notice the higher densities in spiral arms in the inner 
regions in the B ^ 0 case. 


its presence has been linked with a growth in disk eccen- 
tricity (see also [63]). The nonaxisymmetric feature is 
stronger for models with cooling. Hence, a more realistic 
calculation with radiative transfer is necessary to assess 
the strength of nonaxisymmetric structure in a circumbi- 
nary disk. 

As expected, the rest-mass density contrast between 
the two accretion streams becomes larger with smaller 



FIG. 7. Surface density profiles for several q— 1 cases: Ini- 
tial profile (yellow, solid line), B — 0 no-cooling (green circles) 
& cooling (green crosses), B ^ 0 no-cooling (red, solid line) & 
cooling (black, dashed line) cases. Apart from the initial pro- 
file all other profiles are orbit-averaged over the last 10 orbits 
of evolution, beyond which the profile does not change appre- 
ciably. Notice the clear emergence of a pile-up near the inner 
disk edge in the cooling S / 0 case, as well as the relatively 
small change relative to the initial data in the B — 0 cases, 
and the similarity between A = 0 and A ^ 0 unmagnetized 
cases. 


mass ratios. This effect is easily seen when comparing 
the rest mass density contours in the equatorial plane for 
q = 0.25, q = 0.125 and q = 0.1 cases (see Figs. 12 and 
15). 

In the q = 1 and 0.5 no-cooling cases we observe time 
variations in the density of the streams relative to each 
other: for about half an orbit one stream is stronger than 
the other. 

We find that the supply of material channeled onto the 
BHs is sufficient to keep the BHs immersed in a persistent 
gaseous environment with b 2 / po ^ 10 -3 . This means 
that the force-free electrodynamics approximation may 
be inadequate to globally describe the systems considered 
here. 

For q — 0.1 there is hardly a low-density hollow (see 
Fig. 14). This is also revealed by the inner disk edge being 
close to or inside the orbit of the secondary, especially in 
the no-cooling case. In the q = 0 limit no hollow appears, 
as expected. However, we observe a region of lowered 
surface density near and inside the ISCO of the primary 
BH. 


2. Inner disk edge 

In Newtonian 2D studies of geometrically thin disks 
and large binary separations, the location of the inner 
disk edge is found to be roughly twice the binary sepa- 
ration, independent of q ; see, e.g., Table I in [59]. For 
the geometrically thick disks and binary separations we 
are considering, we find the inner disk edge in the re- 
laxed state to be dependent on q and whether cooling is 
enabled. 

In all cases (see the snapshots in Figs. 8-15), the inner 
disk edge is closer than predicted by Newtonian thin- 
disk calculations [59, 63, 71]. In the equal-mass cooling 
case the inner edge remains closer to the initial one (see 
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FIG. 8. Contours at select times of rest-mass density normalized to the initial maximum p max (log scale) in the equatorial 
plane. The plot corresponds to the q— 1 no-cooling case. Here p max — 5.6 x 10“ 11 ( lOM^yf- 1 ) ( 10 8m q ) & cm -3 . 


Fig. 7). The trend is such that the smaller the mass 
ratio the closer the disk edge is to the binary orbit. In 
the q = 0.1 no-cooling case the inner disk edge effectively 
coincides with the orbit of the secondary. We report the 
value for r- m found in each case in Table III, and plot r- in 
vs q in Fig. 17. 


3. Surface density 


In Fig. 18 we show the surface density (E) profiles of 
the relaxed disks, averaged over the last 10 binary orbital 
periods (for all mass ratios, cooling- and no-cooling). For 
all cases the relaxed state deviates strongly from the ini- 
tial profile, highlighting the importance of evolving the 
system for at least a viscous time scale t v i S at the radii 
of interest, where 
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2 R 2 

3 vM 


= 8485 — 
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V0.1 J 


H/R 
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-2 


R 

ISM 


3/2 


( 11 ) 

and where v = 2aP/3po^Kep is the shear viscosity, with 
CKep = (M/R 3 ) 1 / 2 . 


No-cooling: The evolutions for q = 1 and 0.5 are 
similar in terms of their E profiles. The other cases 
(q = 0.25 — 0) yield similar E which extend further in 
than for q = 1 and 0.5. The surface density diagnostic 
clearly demonstrates that the Newtonian thin-disk fea- 
ture that the cavity edge is at twice the binary separation 
and independent of the mass ratio does not hold in this 
class of runs. 

Cooling: For non-zero mass ratios, cooling yields a pile 
up of dense gas near the inner disk edge, which is absent 
in the no-cooling runs and is strongest for the equal-mass 
case. Apart from the pile-up at small radii, all evolutions 
have a rather similar profile at larger radii. In the cooling 
cases the cavity edge is farther out than for no-cooling, 
but still closer than twice the binary separation. The de- 
pendence on mass ratio is weaker than in the no-cooling 
cases, but still in some disagreement with the Newtonian 
thin-disk calculations. 


4- Effective a- stress 

Figure 19 shows the averaged effective Shakura- 
Sunyaev a parameter profiles for all mass ratios, both 
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TABLE III. Table summarizing our main results. Columns show case label, inner disk edge normalized to the binary separation 
rin/d, the mean accretion rate and its approximate standard deviation normalized to the no-cooling single BH mean accretion 
rate ((Mbhbh) =t S Mbhbh) / (M q =o) : main Fourier frequencies normalized by the binary orbital period in the Fourier analysis 
of the accretion rate, Shakura-Sunyaev <a-parameter a , the mean Poynting luminosity (Lem), and the mean cooling luminosity 
(Lcooi) both normalized to the binary mean accretion rate. A dash indicates “no information available”. 


Case 

Tin / d 

((Mbhbh) ± <5Mbhbh)/ (M q =o) 

/ / /orb 

a 

(Lem)/ (Mbhbh)c 2 

(Lcooi)/ (Mbhbh)c 2 

l:lnc-mr 

0.89 

0.43 ± 50% 

(1.0,1.5) 

0.04 

0.013 

- 

l:2nc-lr 

1.33 

0.24 ± 60% 

(0.7, 1.5) 

0.05 

0.012 

- 

l:4nc-lr 

1.24 

0.36 ± 20% 

(0.7, 1.5) 

0.03-0.1 

0.010 

- 

l:8nc-lr 

1.06 

0.41 ± 40% 

0.7 

0.03 - 0.06 

0.010 

- 

l:10nc-lr 

0.92 

0.50 ± 50% 

- 

0.07 

0.017 

- 

Onc-hr 

0.32 b 

1.00 ±24% 

- 

0.05 

0.011 

- 

l:lc-mr 

1.48 

0.43 ± 60% 

(0.5, 1.5) 

0.2 

0.004 

0.127 

l:2c-lr 

1.65 

0.36 ± 30% 

1.0 

0.12 

0.003 

0.110 

l:4c-lr 

1.57 

0.32 ± 60% 

1.0 

0.1 

0.002 

0.107 

l:8c-lr 

1.46 

0.31 ± 30% 

0.6 

0.08 

0.006 

0.096 

l:10c-lr 

1.36 

0.43 ± 30% 

- 

0.013 

0.006 

0.081 

Oc-lr 

0.39 

0.62 ± 10% 

- 

0.012 

0.002 

0.115 


a In some cases we quote a range of values because a single radially averaged value would overestimate a. 

b For ease of comparison, in the single BH case we normalize to 10 M even though it does not correspond to an orbital separation. 


for cooling and no-cooling models in the relaxed state. 
The average is taken over the last 10 binary orbits. In 
Table III we also quote characteristic a values obtained 
by additionally averaging over radii from the inner disk 
edge out to the location of the density maximum. We 
plot a vs q in Fig. 17. 

In all cases we observe larger values for a in the cooling 
cases than for the no-cooling cases, and there is always a 
steep increase in a-stress at smaller radii. For the q — 1 
cooling case we find a(r[ n < r < r max ) ^0.2 and a(r[ n < 
r < r max ) ^ 0.1 for all other cooling cases. A typical 
value for all no-cooling cases is a(r- in < r < r max ) ~ 
0.05. The higher stress for cooling cases results from the 
additional gas compression and associated amplification 
of magnetic fields when cooling is allowed. 

5. Accretion rates 

We compute the accretion rates through the individual 
BHs as well as the total accretion rate onto the binary, 
and normalize these by the (time-averaged) single, non- 
spinning BH accretion rate. 

For a given maximum rest-mass density and total BH 
mass, we find the highest accretion rate in the single BH 
case. This is consistent with the expectation that the 
absence of a tidal-torque barrier will allow matter to flow 
more easily toward the BH(s). 

By contrast, the tidal torque is maximized for q = 1 (all 
else being equal), so the expectation is that the accretion 
rate will be minimum for q = 1. In agreement with this 
expectation we find the accretion rate in the q = 1 cases 
to be smaller than all other mass ratios we consider here. 

In Table III we list the average accretion rate vs mass 


ratio, and plot the results in Fig. 17. The general trend is 
that lower mass ratios have higher accretion rates. In the 
q — 0.1 case the average accretion rate is about 50% that 
of the single BH case with the same initial maximum rest- 
mass density and total BHBH mass, while the average 
accretion rate in the equal-mass case is roughly 33% of 
that in the q = 0 case. 

For q = 1 and 0.5 both black holes accrete at com- 
parable rates whether cooling is applied or not (see four 
upper rows, left panels Fig. 20). However, we observe 
that often the accretion rates on the individual BHs are 
anti-phased, i.e., accretion occurs for half an orbit pri- 
marily on one BH and then for the second half of the 
orbit on the other BH. This behavior in the relaxed state 
is due to an “alternating” pattern in which denser mate- 
rial primarily plunges first through one stream and then 
through the other. 

In Fig. 20 we also plot the accretion rates for q = 0.25 
and 0.1, with and without cooling (three lower rows, left 
panels). It is apparent that for q = 0.25 the accretion 
rate onto the primary is comparable to that onto the 
secondary when A = 0. However, when A % 0 the domi- 
nant contribution to the total accretion rate comes from 
the primary. In the q = 0.1 case we observe that the 
dominant contribution to the total accretion rate comes 
from the primary whether cooling is applied or not, and 
the same holds true for the q = 0.125 case. 

This result in the q — 0.1 case can be qualitatively 
understood through a rough analogy to Bondi-Hoyle- 
Lyttleton accretion: The secondary has a smaller mass 
and moves faster on its orbit reducing its effective cap- 
ture cross section as suggested by Bondi- Hoyle-Lyttleton 
accretion (see e.g. [47]). Also the surface area of the sec- 
ondary is roughly a factor Mg H /m^ h ^ 100 smaller than 


16 


40 

t=5395M 

t=6000M 

t=6720M 

20 





* 




■25 0 25 

X/M 


25 0 25 

X/M 


FIG. 9. Contours at select times of rest-mass density normalized to the initial maximum p m ax (log scale) in the equatorial 


plane. The plot corresponds to the q — 1 cooling case. Here /Wx — 
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FIG. 10. Contours at select times of rest-mass density normalized to the initial maximum p max (log scale) in the equatorial 
plane. The plot corresponds to the q = 0.5 no-cooling case. Here p m ax — 8.3 x 10 -11 ( 15 ^" B y H r i ) ( iq8m 0 ) S cm -3 . 
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FIG. 11. Contours at select times of rest-mass density normalized to the initial maximum p m ax (log scale) in the equatorial 
plane. The plot corresponds to the q = 0.5 cooling case. Here /Wx — 4.2 x 10 -11 ( 2 85M 0 /yr ) ( io^Mq ) S cm -3 ~ 4.2 x 



that of the primary. Note however, that the secondary 
plays a role in stripping matter off the inner disk edge ef- 
fectively as it orbits closest to (or even through) the disk, 
so the accretion rate onto the secondary is not generally 
expected to be 100 times smaller than the accretion rate 
onto the primary. 

In particular, in the q = 0.1 no-cooling case a dense 
persistent structure co-orbits with the secondary (see 
Fig. 14). The density in this structure exceeds the den- 
sity of matter near the primary by more than a factor of 
two. 


6. Variability 

We now report results from the Fourier analysis of the 
accretion rate. These can be seen in the right panels of 
Fig. 20. A summary of the primary Fourier modes for 
the different cases is presented in Table III. 

In the q = 1 case a Fourier analysis of Mbh (accre- 
tion rate onto the primary) and rhbh (accretion rate onto 
the secondary) reveals a characteristic frequency near 
(2/3)M£2bhbh, in agreement with [76]. The analysis of 
A^bhbh (the total binary accretion rate) gives a domi- 
nant Fourier mode with a frequency twice as high. We 
observe a peak at the binary period only for q = 0.5 
and 0.25, and only for the cooling cases (see Fig. 20). 
These results indicate that if the variability in the accre- 
tion rate is directly translated into a variability of EM 
signatures, inferring the binary frequency from EM ob- 
servations may not be straightforward. We observe vari- 
ability for other mass ratios as well. The frequencies for 
the equal mass case also appear in other cases, in addi- 
tion to other weaker contributions, but a clear trend is 
not evident. The most prominent and clean periodic sig- 


nature occurs in the q = 0.5 no-cooling case (see Fig. 20 
and discussion in [64]). Other strong Fourier modes are 
observed in q = 0.5 cooling and the q = 0.25 cases. For 
q = 0.125 and q = 0.1 no significant periodicities are 
observed. 

In the q = 0.1 case the variability is dominated by 
variations in the accretion flow onto the primary. The 
Fourier analysis yields rather irregular accretion, i.e. not 
very pronounced frequencies. The secondary accretes at 
several pronounced frequencies, but the amplitude of the 
variations is much smaller. 

In [64, 65] the dependence on q of the variability 
was studied for geometrically thin (“locally isothermal” 
disks) and proposed as a key feature to observationally 
distinguish accreting BHBHs from standard, single BH 
AGNs. In our Fourier analysis the individual peaks are 
less significant than in [64, 65] and the Fourier spectrum 
yields a more complex structure. This discrepancy is 
likely due to a combination of additional effects including 
differences in the viscosity prescription (i.e. MHD turbu- 
lence vs. a- viscosity ) , cooling prescriptions, thin vs. thick 
disks, 2D vs 3D, and the EOS. These differences result in 
geometrically thin vs thick disks, which are seen to have 
different variability. 


7. Luminosities 

We compute the cooling (L coo i) and Poynting (Lem) 
luminosities, as well as the energy loss rate due to out- 
flowing matter (L gas ) for all cases. We typically find L gas 
is comparable to Lem and quite smaller than L coo i. We 
highlight representative cases and summarize all values 
for the different q characterizing the relaxed state in Ta- 
ble III. 




18 



-25 0 25 

X/M 


25 0 25 

X/M 


FIG. 12. Contours at select times of rest-mass density normalized to the initial maximum p max (log scale) in the equatorial 
plane. The plot corresponds to the q = 0.25 no-cooling case. Here p max — 6.7 x 10 -11 ( inA/'‘"y 1 "- 1 ) ( io^m s ) 6 cm_3 - 
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FIG. 13. Contours at select times of rest-mass density normalized to the initial maximum p max (log scale) in the equatorial 
plane. The plot corresponds to the q = 0.25 cooling case. Here p max — 3.75 x 10 -11 ( 2 27 M 0 /yr ) ( io*m q ) 6 cm -3 ~ 
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FIG. 14. Contours at select times of rest-mass density normalized to the initial maximum p m ax (log scale) in the equatorial 
plane. The plot corresponds to the q — 0.1 no-cooling case. Here p ma x — 6 X 10“ 11 ( hmq H yf- 1 ) ( 10 sm 0 ) S cm -3 . 
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FIG. 15. Contours at select times of rest-mass density normalized to the initial maximum p m ax (log scale) in the equatorial 
plane. The plot corresponds to the q = 0.1 cooling case. Here /Wx — 3.5 x 10 -11 ( 2 85M 0 /yr ) ( io^Mq ) S cm -3 ~ 3.5 x 

10 -11 ( LEdd ) ( lo 8 m q ) § cm -3 . The gas is denser everywhere compared to mass ratios closer to unity; compare to Figs. 8, 


9, 10, 11, 12, 13. 
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FIG. 16. Contours of rest-mass density (linear color scale) normalized to the initial maximum in the equatorial plane at 
t — 10000M. Left panel: q — 0 no-cooling. Right panel: q — 1 no-cooling case. 


The large variability seen in the accretion rate is only 
partially reflected in the cooling luminosities and not re- 
flected in the Poynting luminosities. However, these con- 
clusions need to be confirmed with a self-consistent treat- 
ment involving radiative transfer. 

In all cooling cases we find that the cooling luminosity 
is significantly larger than the Poynting luminosity and 
the energy loss rate due to outflowing matter by almost 
a factor of 10 (see Table III). 

We estimate and compare the contributions to the 
cooling luminosity from various regions in the disk: the 
outer disk, the inner edge, and the cavity (see Fig. 21). 
We find that although the outer disk gives the largest 
contribution, the inner edge and cavity interior are a sub- 
stantial portion (^ 30%) of the total cooling luminosity. 
Therefore, the activity in the cavity cannot be ignored, 


as has been done in earlier studies. 


8. Opacities 

We estimate the Thompson scattering (r es ) and free- 
free absorption (r^) optical depths in all cases. We do 
not find a strong dependence of r* = (SresTff) 1 / 2 , r es 
and Tff on the mass ratio. Our crude analysis shows 
r* - 0{l){L b /L E dd) 9/16 {M/10 8 M Q )- 1 / 16 throughout 
the bulk of the disk. In conjunction with the radia- 
tion pressure dominance found for near Eddington ac- 
cretion rates, this result justifies our choice of T = 4/3 
for the bulk of the disk. Some cases are marginal, how- 
ever given the crudeness of our estimate and scaling ar- 
guments the choice of adiabatic index is adequate. The 
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The lower limit to the equatorial temperature in the 
cavity for all cases (assuming p d e = aT 4 j is 
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FIG. 17. Mean accretion rate M (normalized to the no- 
cooling single BH accretion rate), a, and n n (normalized to 
twice the binary separation) as functions of q for cooling and 
no-cooling cases (all at low resolution). In the absence of 
medium and high resolution runs for g / 1 we place error 
bars in M based on the q — 1 resolution study. These error 
bars are chosen to be 50%, corresponding to the fractional 
difference between the high and low resolutions runs in the 
q = 1 case (see Sec. VA). The error bar in the q = 0 case 
is 15% corresponding to the fractional difference between the 
high and low resolutions runs in the q — 0 case. 


dominant source of opacity in all cases is electron scatter- 
ing. Within the cavity we find that outside the accretion 
streams the matter is optically thin. This means that 
radiation from the cavity can freely stream out, and it 
is likely that (depending on the local temperature) cy- 
clotron lines may give rise to a nonthermal component 
to the emergent EM spectrum. 

We find that r* is affected by cooling. The runs with 
cooling have larger r* than those without cooling. 


9. Characteristic EM radiation frequencies 


The characteristic effective temperatures [see Eq. (9)] 
are 


T eff ~ 10 5 
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The corresponding characteristic thermal radiation fre- 
quencies (z^bb ~ kBT e f[/h) are reduced by a redshift fac- 


implying that the electrons are nonrelativistic, and hence 
they emit cyclotron and not synchrotron radiation. This 
result, too, should be confirmed with radiative transfer 
calculations. Typical cyclotron frequencies in the cavity 
then are 
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Note that due to the large radiation pressure near the 
Eddington limit, one expects that any dust will be blown 
away from the disk and may accumulate at much larger 
radii than our computational domain. This dust is likely 
to absorb the optical/UV radiation and re-emit it in the 
IR [52, 134]. 


1 0. Outflows and jets 

In Fig. 22 (q = 0.1 and 1 cases) we plot the ratio 
6 2 /(2po), which equals the terminal Lorentz factor in ax- 
isymmetric steady-state jet flows. Close to the BH, val- 
ues approaching b 2 / (2 po) ~ 20 are common, dropping to 
b 2 /(2p 0 ) 10 at larger heights Z (in the funnel). We 

observe mildly relativistic outflows in all cases. Our 3D 
visualizations of the B-field lines (see Fig. 23) unambigu- 
ously show that there are field lines emanating from each 
BH horizon and extending into the polar regions. Near 
the BHs [r ~ O(10M)] this leads to a dual jet structure 
that at larger radii [r ^ O(100M)] merges into one com- 
mon helical structure. Due to this effect the dual jets 
may not be detectable individually. In the context of 
force-free simulations around binary black holes, the ex- 
istence of individually detectable dual jets was proposed 
in [80]. However, in [135] it was shown that while a dual- 
jet component is present, it is subdominant with respect 
to the predominantly quadrupolar EM emission, thereby 
casting tight constraints on the detectability of such dual 
jets. Regardless, the “cavity” contains a lot of dense mat- 
ter so that the assumption of the force-free limit of ideal 
MHD may not be applicable. An MHD calculation can 
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FIG. 18. Disk surface density E-profile for different mass ratios as a function of cylindrical radius (normalized to the binary 
separation d ~ 10M). Left panel: no-cooling cases. Right panel: cooling cases. Except for the initial profile all other curves 
correspond to the relaxed state averaged over the last 10 binary orbital periods. In all these B ^ 0 cases the relaxed profile 
deviates strongly from the initial one. In the A = 0 sequence the density reaches further toward small radii with decreasing q. 
In the A/0 sequence, one can see the pile-up near the inner disk edge, which is strongest for large q. 
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FIG. 19. Shakura-Sunyaev a profiles for all mass ratios. Left panel: no-cooling cases. Right panel: cooling cases. Profiles have 
been averaged over the last 10 binary orbital periods. Cooling cases tend to yield larger values. In all cases the stress increases 
inside the cavity. 


resolve this question, and our results suggest that for 
twin jets to be detectable as individual jets, it may re- 
quire either BH spins misaligned with the orbital angular 
momentum or tilted accretion disks [70] . Although for ge- 
ometrically thin disks binary-disk misalignment may be 
unlikely [136], it may be possible for geometrically thick 
disks. 


VI. CONCLUSIONS 

We presented general relativistic magnetohydro dy- 
namic simulations of magnetized circumbinary accretion 


disks onto binary black holes with mass ratios ranging 
from 1:1 to 1:10. We model the disks using a T-law equa- 
tion of state with T = 4/3 - appropriate for optically 
thick, thermal radiation pressure-dominated fluids. We 
focus on a disk near decoupling and perform our compu- 
tations for ~ 10000M (45 binary orbits). We compute 
the disk structure after the disk has reached a quasi- 
relaxed state. This dynamically quasi-steady state is a re- 
sult of binary tidal torques balancing viscous torques aris- 
ing from MHD turbulence triggered by the magnetoro- 
tational instability (MRI). The tidal and viscous torques 
heat the gas, which is expected to radiate and cool. To 
bracket this possibility we perform runs without cooling, 


22 



6.5 7.0 7.5 8.0 



I 


3 



8.5 9.0 9.5 10.0 0 




I 


t/lOOOM 


^fft/ ^ 


BHBH 


FIG. 20. Left panels: Late evolution of mass accretion rate onto the binary (solid, black lines), primary (red, dashed lines), 
and secondary (solid-dotted, yellow lines). Right panels: Fourier transform of the accretion rates shown on the left. From top 
to bottom we show the q — 1 cooling and no-cooling, q — 0.5 cooling and no-cooling, q = 0.25 cooling and no-cooling, and 
q = 0.1 no-cooling cases. 
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FIG. 21. Contribution to L coo \ from the inner cavity and 
outer disk for the q — 1 cooling case. The notation r < Ro 
means the cooling luminosity integral f v Auoy^gd 3 x, where 
V is the volume within a coordinate sphere of radius Ro. Note 
that about 30% of the total cooling luminosity arises within 
Tin, 0 - 


and runs with extremely rapid cooling, adopting an effec- 
tive cooling prescription that resembles a leakage scheme. 

We study how the structure of the accretion flow is al- 
tered for the five different mass ratios. We employ several 
diagnostics, including the accretion rate and its Fourier 
transform, estimates of the electromagnetic luminosity 
and the expected characteristic frequencies of the emer- 
gent electromagnetic radiation, to compare the cases. 

In the equal mass case we find that simulations without 
magnetic fields underestimate the accretion rates, and 
adopting our cooling prescription, the total luminosity of 
the source by two orders of magnitude. This is due to a 
substantial increase of the amount of matter in the ” cav- 
ity” when magnetic fields are present. We also conclude 
that magnetic fields alter the quasi-steady surface den- 
sity profile. Turbulent B-fields lead to more shock heat- 
ing than the binary tidal torques alone, thereby boosting 
cooling luminosities above the values found in unmagne- 
tized disks. 

The surface density profile of the disk is sensitive to 
the mass ratio mainly in the innermost regions. Cooling 
leads to a density enhancement near the inner disk edge 
(a pile-up) which is strongest in the equal mass case. 

We find that for all mass ratios a two-stream accre- 
tion pattern is present. These streams are attached to 
the horizons, with the density close to the horizons being 
among the densest part of the accretion flow. In partic- 
ular, for the q = 0.1 case the material overflows into the 
inner cavity and (partially) refills the hollow present in 
the initial data. During this process a persistent dense 
structure forms around the secondary. This behavior sug- 
gests that it would be inadequate to ignore the flow inside 
the cavity, which only simulations in full GR can treat 


correctly. 

The average binary accretion rates relative to the sin- 
gle, non-spinning BH case range from 50% (q = 0.1) to 
24% (q = 0.5) with a general trend that the accretion 
rate becomes smaller as q — )> 1. For q = 0.5 and q = 1, 
both BHs accrete at a comparable rate (on average), but 
for roughly half a binary period one of the two accretion 
streams is significantly stronger than the other, boosting 
the accretion rate onto one BH and diminishing that onto 
the other. For smaller q the accretion and variability is 
increasingly dominated by the primary BH, especially in 
the cases with cooling. For the single BH (q = 0) signifi- 
cant variability ceases. In general we do not observe ev- 
idence for variability exactly at the binary period. Two 
exceptional cases in which the binary frequency is de- 
tected are the q = 0.5 and 0.25 cooling models. 

We find that the variability in the accretion rate is 
not reflected in the variability of either the Poynting or 
the cooling luminosity. We also find that our cooling 
luminosity is always larger than the Poynting luminos- 
ity, though careful GRMHD simulations with radiative 
transfer are necessary to confirm these findings. 

The cavity radiates an amount comparable to the outer 
disk. Only the innermost regions reflect the strong vari- 
ability seen in M, but in general the radiation from 
the outer disk smears these variabilities out. We ten- 
tatively conclude that it will be challenging to distin- 
guish between single BH and binary BH AGN sources, 
unless other effects such as binary disk misalignment are 
present. However, radiative transfer calculations are re- 
quired to confirm this result. 

We observe nonaxisymmetric structure in the relaxed 
disk. A “lump” forms near the inner disk edge and is 
strongest when we allow for cooling. 

All of our evolutions reveal magnetic field lines ema- 
nating from each of the two horizons, forming dual jets 
which merge into one helical structure above the polar 
regions. 

We estimate that the effective optical depth is 
r* - O(l)(L b /£ Ed d) 9/16 (M/lO 8 M 0 )-Vi6 through- 

out the bulk of the disk, hence the disks we con- 
sider are optically thick to absorption. The char- 
acteristic effective temperature of our disk models 

is T eff ~ lO 5 (L b /L Edd ) 1/4 (M/lO 8 M 0 )“ 1/4 K. Ex- 
pected frequencies of the thermal radiation are v u ~ 

lO 15 (L b /i E dd) 1/4 (M/lO 8 M 0 )- 1/4 /(l + z) Hz (opti- 
cal/near UV). Therefore, instruments such as LSST, 
WFIRST, and PanSTARRS will be ideally suited to 
study these sources. In the cavity, we find that r* < 
1, hence cyclotron lines may be directly observable 
as a nonthermal component of the spectrum. Typi- 
cal cyclotron frequencies in the cavity then are z/ cy ^ 

10 6 (M/lO 8 M 0 )“ 1/2 (L b /L E dd) 1/2 (l + z)- 1 Hz, which 
fall in the radio spectrum. Although considerable cau- 
tion must be applied to these estimates, as we do not 
account for mircophysics and radiative transfer in our 
simulations, they may serve as useful guides to subse- 
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FIG. 22. Contours of b 2 /2po, (log color scale) in a meridional slice, and fluid velocity (arrows) at t ~ 10000M. Left panel: 
q = 1 cooling case. Right panel: q — 0.1 cooling case. 


quent, more detailed investigations and to astronomical 
instruments that may be able to observe the scenarios we 
are simulating. 

In upcoming work we plan to evolve the models pre- 
sented here through inspiral, merger, and the post- 
merger phases. We expect afterglow emission [37, 40] 
similar to [76]. The latter phase will model the process 
by which the disk material viscously diffuses into the cav- 
ity following merger, leading to a brightening of the disk 
[40, 41]. In the future we plan to include more realistic 
cooling and radiative transfer and use higher resolution 
for a more accurate modeling of these sources. It will be 
crucial to determine in this context the level of variability 
in the luminosity, and the EM spectrum of the source. 
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Appendix A: Numerical stability of cooling schemes 
in low density regions 

The evolution equation for the specific thermal energy 
without magnetic fields, but including radiation with an 
effective emissivity A = dU/dr as measured by an ob- 
server comoving with the fluid, is given by [117] 
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^ 6 th _ Pth dpo 1 * / a -i \ 

7 2 j V 

dr pg dr p 0 

where e t h = e — eo, and where e is the total specific 
energy and eo the total specific energy of the fluid ele- 
ment calculated with the EOS at t = 0, i.e., P = KqPq, 
eo = P/po(r — 1). We thus account for cooling by speci- 
fying A. 


1. Cooling prescriptions 

Two effective emissivities that have been adopted in 
the literature are: 


a PoGh 
A i — 5 

t c 

where r c is the cooling time scale, and 


(A2) 
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FIG. 23. Volume rendering of rest-mass density normalized to its initial maximum value (color coding) and magnetic field lines 
for the q = 1 cooling (medium resolution) case. White field lines emanate from the BH apparent horizons. Incipient jets are 
launched from these systems. 
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As we will see below, it is the factor e/eo by which A 2 
differs from Ai that leads to the completely different nu- 
merical stability properties of these two cooling schemes. 

2. Numerical properties of cooling prescriptions 


where S = K = P/po and So is both the initial (un- 
shocked) and target value. We introduced (A2) in [117], 
and (A3) was introduced in [72, 124]. While these two 
schemes achieve similar properties, i.e., drive K back to 
its initial value, they have fundamentally different nu- 
merical stability properties. Here we explain why this is 
so. 

First, let us recast our cooling emissivity Ai in a form 
that is similar to A 2 to show that our cooling prescription 
also drives K to Kq. The pressure is given by 

P = Poe(T - 1) (A4) 

and similarly for the initial pressure 


Po = Po<4)(r - 1). 
The difference between the two is 


(A5) 


A P = p 0 (e - e 0 )(r - 1) = /9oe t h(r - 1) = P 0 — • (A6) 

eo 


But, 


A P_(K- Ko)p$ _AK AS 


= ~T7~ — ~TT- (A7) 


KopI K 0 So 
Combining Eqs. (A6), (A7), and (A2) yields 
Poe o AS 


Ai = 


r c So 


(A8) 


Insert our cooling prescription in Eq. (AI) and drop 
the first term on the right-hand-side (RHS) of Eq. (AI), 
i.e, assume that no adiabatic compression or expansion 
takes place. Then we obtain 


dtth ^th 

dr r c ’ 


(All) 


i.e., exponential cooling of the excess thermal energy with 
cooling time scale r c . For an explicit numerical scheme 
to be Courant stable we must set the maximum timestep 
At < r c . To see this we can we use a simple Euler explicit 
integration scheme to write Eq. (All) in finite difference 
form as 


f n+1 - e 
6 th 6 


th 


Hh 


At t c 

Then the amplification factor is given by 


f = 


Tn+1 

e th 


c th 


1 


At 


(A12) 


(A13) 


For numerical stability the magnitude of the amplifica- 
tion factor must be less than 1 


i/i < i 


i - 


At 


< 1 =>• At < 2 t c . 


(A14) 


By contrast, inserting A 2 [Eq. (A10)] in Eq. (AI) and 
dropping the first term on the RHS of Eq. (AI) yields 


de th _ _ 2 fth _e_ 
dr t c e 0 ’ 


(A15) 


The emissivity A 2 for AS > 0, i.e., the only case when 
the emissivity is not 0, becomes 


A 2 


2 poe AS 
t c So 


(A9) 


Thus, the two prescriptions Eqs. (A8) and (A9) are 
similar and both drive the gas to constant initial entropy 
So- But the two prescriptions are not the same. Another 
way to see this is to write A 2 in a form that resembles 
Ai. Using Eqs. (A6), (A7), we can write (A9) as 


A 2 


2ppeth e 
Tc to' 


(A10) 


hence the effective cooling time scale is now 

Teff = \ ~e Tc ' *- A16 ) 

We now require At < r e ff for numerical stability. How- 
ever, low density regions can become strongly shock 
heated, in which case eo « e, implying r e ff r c . As 
a result, a violation of the Courant stability condition 
can quickly occur that will terminate the computations, 
if At < t c is used as the stability condition. This is 
a disadvantage of the A 2 prescription, as it means that 
one must impose a density cutoff for the simulations to 
be stable and retain relatively large timesteps. This is 
something we do not need to require with our cooling 
emissivity Ai, as the fixed cooling time scale r c is, at the 
same time, the effective cooling time scale. The difference 
is important as we want to probe low density regimes in 
the disk without artificial cutoffs. 
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